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Abstract 

Insect population decline has been reported worldwide, including those of pollinators important for eco- 
system services. Therefore, conservation actions which rely on available rigorous species distribution data 
are necessary to protect biodiversity. Niche modeling is an appropriate approach to distribution maps, but 
when it comes to bumble bees, few studies have been performed in South America. We modeled ecological 
niches of nine Colombian Bombus species with MAXENT 3.4 software using bioclimatic variables available 
from WorldClim. This resulted in maps for each species that show the potential distribution area at the pre- 
sent time. Modeled species maps accurately represent potential niches according to the description of bio- 
climatic conditions in the species’ habitat. We grouped the species into three clusters based on our results, 
as well as on distributional information from literature on the topic: High Mountain, Mid- Mountain and 
inter-Andean, and the Amazon and Eastern Plains Basin. Niche modeling depicted bumble bee species’ dis- 


tribution in Colombia, the results of which can serve as a useful tool for conservation policies in the country. 
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Introduction 


Bumble bee species (Bombus spp.) are one of the widest studied taxa of bees in the 
world and Colombia regarding their ecological traits, the genetic composition of 
populations, pathogens, and their role in pollination services (Cure and Rodriguez 
2007; Gamboa et al. 2015; Jaramillo-Silva et al. 2018; Lotta-Arevalo et al. 2020). Buzz 
pollination carried out by bumble bees allows for the pollination of several crops of 
economic importance in Colombia, such as Solanaceae like potatoes, tomatoes, pa- 
prika, etc. (Cure and Rodriguez 2007; Riafo et al. 2015; Vinicius-Silva et al. 2017) or 
Passifloraceae like passion fruit, gulupa, curuba, etc. (Pinilla-Gallego et al. 2016). Also, 
bumble bees are important for the pollination of many plant species in natural ecosys- 
tems like paramos and cloud forests, which are key for water regulation in Colombia 
(Rubio 2012). For example, they serve in the pollination of several species of “frailejon” 
of the genus Espeletia (Fagua and Bonilla Gémez 2005). 

However, a drastic decline in insects has been documented worldwide (Hallmann 
et al. 2017; Seibold et al. 2019; Klink et al. 2020). A massive decline in the species dis- 
tribution of most bumble bees has been reported, especially in developed regions like 
Europe and North America (Goulson et al. 2008). In Europe, there have been reports 
of declines since the 50’s. Just in the United Kingdom, three species out of twenty-five 
have gone extinct, and eight report major declines (Goulson 2003). Similar results 
have been reported in North America. In Illinois alone, four species have become lo- 
cally extinct (Grixti et al. 2009). There is no population trend research of this kind 
performed in Colombia, and as one of the richest countries in terms of biodiversity 
(Rangel 2015), the country should increase research on its pollinators in order to de- 
velop preventive management of natural resources. 

Biogeographic information about species is vital for conservation planning. The 
study of species distribution patterns has two complementary approaches: historical 
biogeography, which elucidates the causal mechanisms of current distribution, and 
ecological biogeography, which evaluates the ecogeographic factors that are currently 
shaping the distribution of species, looking for patterns in characteristics that are re- 
quired for a species’ long-term survival. It looks at abiotic characteristics like humid- 
ity, temperature, or salinity and ecological characteristics like interactions with other 
organisms or genetical features. Research in this last approach is usually performed at 
a local scale a short time frame (Pérez-Malvaez and Gutiérrez 2003). 

Niche modeling is one of the methods used in ecological biogeography. It repre- 
sents the fundamental niche without considering the realized niche, which requires 
detailed field research and confirmed species samples. Niche modeling aims to find the 
potential distribution area of a species (Soberon and Miller 2009). It is important to 
have good quality and quantity of information for the models, as they can be biased if 
the data for a species is scarce. 

Thus, niche modeling can provide potential distribution maps for bumble bees in 
Colombia, updating and complementing the previous available maps, made by Abram- 
ovich et al. (2004). This information could be an important tool to use in both future 
studies with an evolutionary and ecological aim, as well as in conservation efforts by 


Niche modeling of Bumblebee species in Colombia p56} 


informing environmental conservation plans, political decisions, and public policies in 
Colombia, and by strengthening existing efforts like the Colombian Pollinators Initiative 
(CPI) (Nates-Parra 2016). As the importance of bee conservation continues to escalate, 
conservation actions are necessary to protect bees’ biodiversity (Sarospataki et al. 2005). 


Materials and methods 


Study area 


The geographical extent used for each species in the modeling process covers the entire 
continental Colombian territory, from the north (13°23.73'N) to the south (4°13.75'N) 
and from the west (81°44.13'W) to the east (66°50.63'W). No buffer was used. 

There are nine species from the genus Bombus reported for Colombia: B. pauloensis 
Friese, 1913 (formerly B. atratus Franklin, 1913), B. excellens Smith, 1879, B. funebris 
Smith, 1854, B. hortulanus Friese, 1904, B. melaleucus Handlirschi, 1888, B. pullatus 
Franklin, 1913, B. robustus Smith, 1854, B. rubicundus Smith, 1854, and B. transversalis 
Oliver, 1789. Across the paramo ecosystem, four species can be found with different 
altitudinal distributions: B. funebris between 2500 and 4750 m, B. hortulanus between 
2100 and 3600 m, B. robustus between 2100 and 3800 m, and B. rubicundus between 
2500 and 3900 m (Pinilla-Gallego et al. 2016). Along the sub-Andean forest or the 
cloud forest, two species can be found: B. excellens between 1500 and 2600 m and 
B. melaleucus between 450 and 2100 m. The most abundant species in the low moun- 
tain strata is B. pauloensis. However, it has a wide altitudinal range between 150 and 
3600 m (Liévano et al. 1991). Finally, two different species can be found across tropical 
forests with warm and wet climates: B. pullatus between 120 and 3500 m (especially in 
the foothills) and B. transversalis between 180 and 1100 m (restricted to the amazon 
forest and its foothills). 


Occurrence data 


We obtained occurrence points from the Wild Bee Research Lab of the Universidad Na- 
cional de Colombia, Bogota (LABUN), the Insect Collection of Universidad del Quindio 
(CIUQ), the Francisco Luis Gallego Entomological Museum (MEFLG), Universidad 
Nacional de Colombia, Medellin, the Bee Collection of the Universidad Militar Nueva 
Granada, the Entomological Museum of the School of Agronomy at Universidad Na- 
cional de Colombia (UNAB), and the database from the Global Biodiversity Information 
Facility (GBIF). Occurrence points were adjusted using the altitude reference reported 
by Pinilla-Gallego et al. (2016). All occurrence point data were originally collected from 
1938 until 2020, throughout different seasons and derived from various methods (See 
Suppl. material 1). The senior author checked all occurrence records to prevent recogniz- 
able errors in georeferencing and taxonomy. The number of points used for each specie 
was: B. excellens 46; B. funebris 150; B. hortulanus 1130; B. melaleucus 47; B. pauloensis 
2128; B. pullatus 549; B. robustus 243; B. rubicundus 426; B. transversalis 43 (Fig. 1). 
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Figure |. Distribution of the points used in the modeling of each specie. 
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Environmental variables 


We used 19 environmental data layers for modeling (Table 1), which were downloaded 
from the WorldClim V. 2.1 website (Fick and Hijmans 2017). The layers represent the 
climate average from the year 1970 to the year 2000 at a spatial resolution of 30 Arc 
seconds (0.93 km”). 


Niche models 


QGIS 2.8 (QGIS Development Team 2016) was used to clip and convert environ- 
mental layers, enabling their use in MAXENT. Niche models were developed using 
MAXENT 3.4 (Phillips et al. 2020). Specific settings were set at the following defaults: 
a maximum of 500 iterations, 10% test point (CrossValidate), without extrapolation 
and cumulative, keeping a limit convergence of 0.00001 and prevalence of 0.5, maxi- 
mum number of background points 10000. 

Two different statistical analyses were performed per species: the Jackknife test to 
evaluate the weight or importance of each variable (Timana de la Flor and Romero 2015; 
Phillips and AT&T Research 2017), and the AUC (Area Under the Curve) test. Jack- 
knife test is a method for validating the samples and model. It shows the representative 
variables for modeling each specie (Shcheglovitova and Anderson 2013). According to 
the results of the Jackknife test, we selected the three variables with the highest percent- 
age of importance, whose values depend on each species. To obtain the AUC, the specific 
settings given above were chosen according to established knowledge about bumble bees’ 
altitudinal distribution (Pinilla-Gallego et al. 2016). Subsequently, an ROC (Receiver 
Operating Characteristics) curve was created, and the AUC was calculated. The AUC 
value reflects the model’s accuracy or its capacity for prediction. The AUC value moves 
between 0 and 1, with 1 representing a perfect prediction. Values over 0,9 are considered 
strong (Pliscoff and Fuentes-Castillo 2011). Finally, ArcMap 10.5 (ESRI 2015) was used 


for reclassifying and visually representing the ecological niche modeling results. 


Table |. Bioclimatic variables included in modeling were obtained from WorldClim. 


Temperature Precipitation 
1 Annual Mean Temperature 12 Annual Precipitation 
2 Mean Diurnal Range (Mean of monthly (max. temp.-min. temp.)) 13 Precipitation of Wettest Month 
3 Isothermality (BIO2/BIO7) (x100) 14 Precipitation of Driest Month 
4 Temperature Seasonality (standard deviation x 100) 15 Precipitation Seasonality (Coefficient of Variation) 
5 Max. temperature of Warmest Month 16 Precipitation of Wettest Quarter 
6 Min. Temperature of Coldest Month 17 Precipitation of Driest Quarter 
7 Temperature Annual Range (5-6) 18 Precipitation of Warmest Quarter 
8 Mean Temperature of Wettest Quarter 19 Precipitation of Coldest Quarter 


9 Mean Temperature of Driest Quarter 
10 Mean Temperature of Warmest Quarter 
11 Mean Temperature of Coldest Quarter 
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Results 


The models showed good results for all nine bumble bee species, with AUC values above 
0.9, although the number of occurrences for Bombus excellens, Bombus melaleucus, and 
Bombus transversalis was low. Bombus pauloensis and Bombus hortulanus had the highest 
number of occurrences, being common in collections (Table 2). Based on our results 
and distributional information from scientific literature, we divided the bumble bee 
species into three distributional groups: high mountain bumble bees, mid-mountain/ 
inter-andean bumble bees, and Amazon and eastern Plains basin bumble bees. 


High Mountain bumble bees 


The potential distribution areas for B. funebris, B. hortulanus, B. robustus, and 
B. rubicundus occurred only in high mountain departments, along a range of different 
small and fragmented areas, with a high probability of occurrence in the central part 
of the Eastern Andes Range (Fig. 2). B. funebris represented the most restricted and 
fragmented species. B. hortulanus occurrences had the lowest altitudinal points for this 
group. B. robustus data was similar to B. hortulanus but included some areas in the 
northern part of the Eastern and the Central Ranges. B. rubicundus showed similar 
potential distribution to B. robustus but with more restricted and fragmented areas. 


Mid-Mountain and Inter-Andean bumble bees 


This group of species showed a wider altitudinal and potential distribution than 
the high mountain species, along and between the three mountain ranges (Fig. 3). 
B. excellens showed a higher probability along most of the Central and Eastern moun- 
tain ranges, as well as in the mid-mountain area of the Sierra Nevada de Santa Marta. 
For B. melaleucus, its potential distribution was indicated along the three moun- 
tain ranges. B. pauloensis exhibited a continuous distribution over the Andean region. 
B. pullatus was concentrated along low altitudinal areas of the inter-Andean valleys. 


Table 2. Bumble bee species data for the area under the curve (AUC obtained by niche modeling), num- 
ber of occurrences used for modeling, and the distributional group selected for each species according to 
our results and previous results by Pinilla-Gallego et al. (2016). 


Species AUC 1970-2000 Distributional group 
Bombus excellens 0.951 Mid-mountain and inter-Andean 
Bombus funebris 0.976 High Mountain 
Bombus hortulanus 0.977 High Mountain 
Bombus melaleucus 0.965 Mid-mountain and Inter-Andean 
Bombus pauloensis 0.911 Mid-mountain and Inter-Andean 
Bombus pullatus 0.917 Mid-mountain and Inter-Andean 
Bombus robustus 0.983 High Mountain 
Bombus rubicundus 0.973 High Mountain 


Bombus transversalis 0.999 Amazon and Eastern Plains Basin 
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Figure 2. Potential distribution maps for High Mountain bumble bees A B. funebris B B. hortulanus 
C B. robustus D B. rubicundus. 
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Figure 3. Potential distribution maps for Mid-Mountain and Inter-Andean bumble bees. A B. excellens 
B B melaleucus C B. pauloensis D B. pullatus. 
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Amazon and Eastern Plains Basin bumble bees 


Bombus transversalis is the only species with a potential distribution along the Ama- 
zon Forest, the Andean foothills, and the Orinoquia region. The eastern part of Meta 
presented the highest altitudinal points for the species at 1150 m. The northern area 
of the biogeographic region of Chocé showed optimal bioclimatic conditions for 
B. transversalis, but no occurrence point was found or used there for this model (Fig. 4). 
The variables with the highest weight for each species are presented in Table 3. 
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Figure 4. Potential distribution map for Amazon and Eastern Plains Basin bumble bees, B. transversalis. 
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Table 3. Bioclimatic variables had the highest weight in the model for each species, according to the 
results obtained by the Jackknife test. 


Species (Distributional group) Variables with the highest weight (Percent of contribution) 
1 7nd 3M 
Bombus excellens (Mid- 8. Mean Temperature of 1. Annual Mean Temperature 6. Min. Temperature of Coldest 
Mountain and inter-Andean) Wettest Quarter 94.35% 94.27% Month 94.04% 
Bombus funebris 4. Temperature Seasonality 6. Min. Temperature of Coldest 1. Annual Mean Temperature 
(High Mountain) 95.50% Month 94.54% 94.03% 
Bombus hortulanus 11. Mean Temperature of 1. Annual Mean Temperature 10. Mean Temperature of 
(High Mountain) Coldest Quarter 96.72% 96.63% Warmest Quarter 96.61% 
Bombus melaleucus (Mid- 5. Max. Temperature of 11. Mean Temperature of 1. Annual Mean Temperature 
Mountain and inter-Andean) Warmest Month 95.20% Coldest Quarter 95.16% 95.14% 
Bombus pauloensis 6. Min. Temperature of Coldest | 11. Mean Temperature of 8. Mean Temperature of 
(High Mountain) Month 91.30% Coldest Quarter 91.18% Wettest Quarter 91.09% 
Bombus pullatus (Mid- 11. Mean Temperature of 8. Mean Temperature of 4. Temperature Seasonality 
Mountain and inter-Andean) Coldest Quarter 83.91% Wettest Quarter 83.53% 76.67% 
Bombus robustus 5. Max. Temperature of 6. Min. Temperature of Coldest 8. Mean Temperature of 
(High Mountain) Warmest Month 93.74% Month 93.72% Wettest Quarter 93.51% 
Bombus rubicundus 6. Min. Temperature of Coldest 1. Annual Mean Temperature 8. Mean Temperature of 
(High Mountain) Month 96.70% 96.70% Wettest Quarter 96.48% 
Bombus transversalis (Amazon 13. Precipitation of Wettest 11. Mean Temperature of 6. Min. Temperature of Coldest 
and Eastern Plains Basin) Month 98.88% Coldest Quarter 97.87% Month 96.94% 
Discussion 


These potential distribution maps for bumble bees in Colombia improve the previous 
maps available (Abrahamovich et al. 2004) for the species, with more precise areas of 
potential distribution and a visual schema for patterns already described in the litera- 
ture (Pinilla-Gallego et al. 2016). 

The four high-mountain species are associated with the paramo and high Ande- 
an ecosystems (Pinilla-Gallego et al. 2016). Paramo ecosystems are characterized as 
biogeographical islands, highly isolated areas with a great deal of endemism and low 
genetic flow between populations (Lotta-Arevalo et al. 2020). As a result, the poten- 
tial distribution of these species is highly fragmented. B. excellens, B. melaleucus, and 
B. pullatus are associated with mountain forests, and it is remarkable that, even though 
they showed a wide suitable area (Fig. 3A, B, D, respectively), there is a low number of 
occurrences. This indicates a small population size and high vulnerability. For example, 
B. excellens can only be found in the cloud forest, an ecosystem with an accelerated rate 
of deforestation (Nates-Parra 2006). B. pauloensis (Fig. 3C) shows a wide altitudinal 
and potential distribution consistent with its high plasticity in the habitat chosen for 
nesting (Liévano et al. 1991; Nates-Parra et al. 2006), and this is reflected by its wide 
distribution along the mountain areas. Its potential distribution will help conservation 
planning, as B. pauloensis positively impacts national agricultural productivity by pol- 
linating staple foods such as fruits and vegetables (Cure and Rodriguez 2007; Riafo et 
al. 2015; Poveda et al. 2018) due to its adaptability to disturbed habitats. According 
to its biogeographical history, for B. transversalis the Andes Mountain Range repre- 
sents a physical barrier for emigration. Thus, the potential distribution area in the 
biogeographical region of Chocé is not likely (Abrahamovich et al. 2004). The eastern 
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area of Meta stands out as an expansion of previously recorded areas to include the low 
areas of the Andean foothills. 

The most important bioclimatic variables for bumble bees in Colombia are related 
to temperature (Table 3). Thus, drastic temperature changes put bumble bees in a 
condition vulnerable to threat under a climate change scenario (Gonzalez et al. 2021). 
Also, the potential distribution for most species overlaps with densely populated areas. 
Therefore, these maps are a tool for the detailed location of bumblebee biodiversity for 
conservation planning as bumble bee’ populations and their associated services will 


probably soon be reduced (Williams and Osborne 2009; Pinilla-Gallego et al. 2016). 


Conclusion 


The potential distribution for Colombian bumble bee species is reported here. Seven of 
them show a restricted distribution, shaped mainly by temperature restrictions. These 
results are a direct contribution to knowledge about Colombian bumble bee species, 
constituting useful knowledge for conservation, territorial planning, protection plans, 
and environmental management. Thus, by obtaining the most suitable areas for a spe- 
cies, this study can provide the ideal location to breed and reproduce the native species. 
These species contribute to improving the agricultural productivity of several regions 
by pollination. Likewise, this information can help avoid the introduction of foreign 
species for this purpose (Nates-Parra 2016). It is necessary to increase research efforts 
on bees for them to be included in conservation planning. 
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